Diagonally Implicit Runge-Kutta Research Articles

A high-order finite difference method for moving immersed domain boundaries and material interfaces

We present a high-order sharp treatment of immersed moving domain boundaries and material interfaces, and apply it to the advection-diffusion equation in two and three dimensions. The spatial discretization combines dimension-split finite difference schemes with an immersed boundary treatment based on a weighted least-squares reconstruction of the solution, providing stable discretizations with up to sixth order accuracy for diffusion terms and third order accuracy for advection terms. The temporal discretization relies on a novel strategy for maintaining high-order temporal accuracy in problems with moving boundaries that minimizes implementation complexity and allows arbitrary explicit or diagonally-implicit Runge-Kutta schemes. The approach is broadly compatible with popular PDE-specialized Runge-Kutta time integrators, including low-storage, strong stability preserving, and diagonally implicit schemes. Through numerical experiments we demonstrate that the full discretization maintains high-order spatial and temporal accuracy in the presence of complex 3D geometries and for a range of boundary conditions, including Dirichlet, Neumann, and flux conditions with large jumps in coefficients.

Construction of A-stable explicit last-stage diagonal implicit Runge–Kutta (ELDIRK) methods

AbstractELDIRK methods are defined to have an Explicit Last stage in the general Butcher array of Diagonal Implicit Runge-Kutta methods, with the consequence, that no additional system of equations must be solved, compared to the embedded RK method. Two general formulations for second- and third-order ELDIRK methods have been obtained recently in Mahnken [21] with specific schemes, e.g. for the embedded implicit Euler method, the embedded trapezoidal-rule and the embedded Ellsiepen method. In the first part of this paper, we investigate some general stability characteristics of ELDIRK methods, and it will be shown that the above specific RK schemes are not A-stable. Therefore, in the second part, the above-mentioned general formulations are used for further stability investigations, with the aim to construct new second- and third-order ELDIRK methods which simultaneously are A-stable. Two numerical examples are concerned with the curing for a thermosetting material and phase-field RVE modeling for crystallinity and orientation. The numerical results confirm the theoretical results on convergence order and stability.

Diagonally implicit Runge–Kutta schemes: Discrete energy-balance laws and compactness properties

Abstract We study diagonally implicit Runge–Kutta (DIRK) schemes when applied to abstract evolution problems that fit into the Gelfand-triple framework. We introduce novel stability notions that are well-suited to this setting and provide simple, necessary and sufficient, conditions to verify that a DIRK scheme is stable in our sense and in Bochner-type norms. We use several popular DIRK schemes in order to illustrate cases that satisfy the required structural stability properties and cases that do not. In addition, under some mild structural conditions on the problem we can guarantee compactness of families of discrete solutions with respect to time discretization.

Entropy dissipative higher order accurate positivity preserving time-implicit discretizations for nonlinear degenerate parabolic equations

We develop entropy dissipative higher order accurate local discontinuous Galerkin (LDG) discretizations coupled with Diagonally Implicit Runge–Kutta (DIRK) methods for nonlinear degenerate parabolic equations with a gradient flow structure. Using the simple alternating numerical flux, we construct DIRK-LDG discretizations that combine the advantages of higher order accuracy, entropy dissipation and proper long-time behavior. We theoretically prove the entropy dissipation of the implicit Euler-LDG discretization without any time-step restrictions when no positivity constraint is imposed. Next, in order to ensure the positivity of the numerical solution, we use the Karush–Kuhn–Tucker (KKT) limiter, which achieves a positive solution by coupling the positivity preserving KKT conditions with higher order accurate DIRK-LDG discretizations using Lagrange multipliers. In addition, mass conservation of the positivity-limited solution is ensured by imposing a mass conservation equality constraint to the KKT equations. Under a time step restriction, the unique solvability and entropy dissipation for implicit first order accurate in time, but higher order accurate in space, positivity-preserving LDG discretizations with periodic boundary conditions are proved, which provide a first theoretical analysis of the KKT limiter. Finally, numerical results demonstrate the higher order accuracy and entropy dissipation of the positivity-preserving DIRK-LDG discretizations for problems requiring a positivity limiter. In addition, we can observe from the numerical results that the implicit time-discrete methods alleviate the time-step restrictions needed for the stability of the numerical discretizations, which improves computational efficiency.

On the accuracy, stability and computational efficiency of explicit last-stage diagonally implicit Runge–Kutta methods (ELDIRK) for the adaptive solution of phase-field problems

Runge–Kutta algorithms with adaptive step size control provide reliable tools for the solution of initial value problems. Most widespread within the vast family of Runge–Kutta methods are the diagonally implicit Runge–Kutta (DIRK) methods. In particular, the new low-order explicit last-stage diagonally implicit Runge–Kutta (ELDIRK) methods are investigated, combining implicit schemes with an additional explicit evaluation as an explicit last stage. This results in Butcher tableaus with two solutions of different convergence orders suitable for embedded methods, where the higher-order solution is achieved by additional explicit evaluations. Thus, the iterative solution of non-linear systems is omitted for the additional stage, presenting a major reduction in computational cost for the determination of a local error estimate. The key contribution is the application of the novel Butcher tableaus to phase-field problems, solved with the finite-element method, leading to substantial numerical investigations with an efficient approach for diagonally implicit Runge–Kutta schemes. The most important aspects are the extension towards fourth-order methods, the study of the convergence orders for the ELDIRK schemes, their respective stability regions and computational efficiency. A local error estimator is presented capturing the evolution of the phase-field problems, such that adaptive step size control for the new low-order embedded schemes based on an empirical approach for error estimation is achieved. A suitable parallel algorithm is presented with conclusive two-dimensional phase-field simulations based on a Kobayashi–Warren–Carter model including benchmarks for computational efficiency. The higher-order convergence suggested by the novel schemes is confirmed, and their effective results are demonstrated, resulting in a valuable semi-explicit addition to the family of Runge–Kutta time integration schemes.

Numerical investigations of new low‐order explicit last stage diagonal implicit Runge–Kutta schemes with the finite‐element method

AbstractInitial value problems can be solved efficiently by means of Runge–Kutta algorithms with adaptive step size control. Diagonally implicit Runge–Kutta (DIRK) methods are the most popular class among the diverse family of Runge–Kutta algorithms. In this paper, the novel class of low‐order explicit last‐stage diagonally implicit Runge–Kutta (ELDIRK) methods are explored, which combine implicit schemes with an additional explicit evaluation as an explicit last stage. ELDIRK Butcher tableaus are used to control embedded RK methods to obtain solutions of different orders. The lower‐order solution is obtained by classical implicit RK stages and the higher‐order solution is obtained by additional explicit evaluation. As a result, a significant reduction in computational cost is achieved by skipping the iterative solution of nonlinear systems for the additional step. The examination of the heat problem and the use of the innovative Butcher tableau in the finite‐element method are the main contributions of this work. Thus, it is possible to establish adaptive step size control for the new low‐order embedded methods based on an empirical method for error estimation. Two‐dimensional simulations are used to show an appropriate algorithm for the ELDIRK schemes. The new Runge–Kutta schemes' predictions of higher‐order convergence are confirmed, and their successful outcomes are illustrated.

Design of DIRK schemes with high weak stage order

Runge-Kutta (RK) methods may exhibit order reduction when applied to certain stiff problems. While fully implicit RK schemes exist that avoid order reduction via high-stage order, DIRK (diagonally implicit Runge-Kutta) schemes are practically important due to their structural simplicity; however, these cannot possess high stage order. The concept of weak stage order (WSO) can also overcome order reduction, and it is compatible with the DIRK structure. DIRK schemes of WSO up to 3 have been proposed in the past, however, based on a simplified framework that cannot be extended beyond WSO 3. In this work a general theory of WSO is employed to overcome the prior WSO barrier and to construct practically useful high-order DIRK schemes with WSO 4 and above. The resulting DIRK schemes are stiffly accurate, L-stable, have optimized error coefficients, and are demonstrated to perform well on a portfolio of relevant ODE and PDE test problems.

Accuracy and Stability Analysis of the Semi-Lagrangian Method for Stiff Hyperbolic Relaxation Systems and Kinetic BGK Model

In this paper, we develop a family of third order asymptotic-preserving (AP) and asymptotically accurate (AA) diagonally implicit Runge–Kutta (DIRK) time discretization methods for the stiff hyperbolic relaxation systems and kinetic Bhatnagar–Gross–Krook (BGK) model in the semi-Lagrangian (SL) setting. The methods are constructed based on an accuracy analysis of the SL scheme for stiff hyperbolic relaxation systems and kinetic BGK model in the limiting fluid regime when the Knudsen number approaches 0. An extra order condition for the asymptotic third order accuracy in the limiting regime is derived. Linear von Neumann stability analysis of the proposed third order DIRK methods are performed to a simplified two-velocity linear kinetic model. Extensive numerical tests are presented to demonstrate the AA, AP, and stability properties of our proposed schemes.

Implementation of Artificial Compressibility Method for Steady and Unsteady Incompressible Flows Based on Finite Volume Method of Unstructured Grids

The development of a new computational fluid dynamics solver for simulating incompressible flows is described. Unlike traditional pressure-based solvers, the artificial compressibility method simulates both steady and unsteady states. Some classical numerical examples are used to validate the solvers. These findings show that the artificial compressibility method’s LES/DNS and steady flow simulations are competitive with the best previously reported methods. The following characteristics make this paper unique: (1) We implement the artificial compression method in OpenFOAM for the first time and show the detailed mathematical process of the algorithm so that readers can easily reproduce and improve it; (2) steady flow simulation case, the convergence time of the present method is 52% of the pressure base solver; (3) to accelerate the convergence rate of each time step to the incompressible state, the L-stable Singly Diagonal Implicit Runge-Kutta method is introduced into the pseudo-time advance, which enlarges the limit of maximum Courant number to accelerate convergence.

Efficiency of diagonally implicit Runge-Kutta time integration schemes in incompressible two-phase flow simulations

Following a rapidly growing popularity of the use of computational fluid dynamics (CFD) for free surface flow applications, especially in the naval field, cost efficient numerical schemes are highly desirable. In the present study, diagonally implicit Runge-Kutta (DIRK) time integration schemes up to third-order are implemented in OpenFOAM for incompressible two-phase flow simulations. The coupling of these time integration schemes with solving the incompressible Navier-Stokes equations and the Volume of Fluid (VOF) convection equation is described in detail. After detailed convergence study, the efficiency of these schemes is evaluated in a number of one-phase and two-phase benchmark cases in terms of their accuracy and CPU-to-accuracy ratio, especially in ranges of space and time resolution parameters of practical use. It is shown in particular that the proposed second order Runge-Kutta schemes present the best CPU-to-accuracy compromise among the tested schemes.

Stability Analysis of Diagonally Implicit Two Derivative Runge-Kutta methods for Solving Delay Differential Equations

The stability properties of fourth and fifth-order Diagonally Implicit Two Derivative Runge-Kutta method (DITDRK) combined with Lagrange interpolation when applied to the linear Delay Differential Equations (DDEs) are investigated. This type of stability is known as P-stability and Q-stability. Their stability regions for (λ,μ∈R) and (μ∈C,λ=0) are determined. The superiority of the DITDRK methods over other same order existing Diagonally Implicit Runge-Kutta (DIRK) methods when solving DDEs problems are clearly demonstrated by plotting the efficiency curves of the log of both maximum errors versus function evaluations and the CPU time taken to do the integration.

A high-order unified stochastic particle method based on the Bhatnagar-Gross-Krook model for multi-scale gas flows

Stochastic particle methods, such as Direct Simulation Monte Carlo (DSMC), are successfully used to simulate rarefied gas flows. However, since traditional stochastic particle methods employ the operator splitting scheme, i.e., decoupling the molecular motion and collisions, they are of low order accuracy in the continuum regime. In contrast, the discrete velocity method can easily construct high order schemes, but its computational and memory cost is quite high due to the huge number of discrete velocity nodes, especially for the high Mach number flows. To improve the efficiency and accuracy of the simulation for multi-scale flows, a high order stochastic particle BGK method was developed in this paper. The DIRK (Diagonally Implicit Runge–Kutta) and WLS-ENO (weighted-least-squares based essentially non-oscillatory) schemes were used to achieve a high order of accuracy in time evolution and space interpolation, respectively. However, the implementation of DIRK and WLS-ENO with computational particles is not straightforward due to the statistical nature of the stochastic particle methods. To circumvent these difficulties, firstly, the micro-macro decomposition of the relaxation term of the BGK model, which was proposed in the unified stochastic particle BGK (USP-BGK) method, was employed. The fluid limit part of the relaxation term was treated following the DIRK scheme, while the kinetic part was solved with an exact integration scheme similar to the traditional stochastic particle BGK (SP-BGK) method. Therefore, in the continuum regime, the present stochastic particle method has the same asymptotic preserving and asymptotic accurate properties as the high order semi-Lagrangian method; and in the rarefied regime, it reduces to the traditional SP-BGK method. Additionally, to overcome the effect of thermal noise in the space interpolation, a new discontinuity indicator, which considers the influence of the thermal fluctuation, was also proposed for the WLS-ENO scheme. Several 1D and 2D benchmark problems with different Mach and Knudsen numbers were tested and analyzed, such as shock and sine wave interaction, Taylor-Green vortex flow and 2D Riemann problem. Compared to the traditional SP-BGK method, this proposed method can achieve a higher order of accuracy and significantly improve the computational efficiency in the continuum regime.

High-order implicit weighted compact nonlinear scheme for nonlinear coupled viscous Burgers’ equations

In this paper a high-order implicit weighted compact nonlinear scheme for nonlinear coupled viscous Burgers’ equations is presented. The fifth-order weighted compact nonlinear scheme is used for the spatial discretization, while the third-order diagonal implicit Runge–Kutta method is used for the time discretization. The generated nonlinear system is solved by the Jacobian-free Newton–Krylov nonlinear solver, which is composed of the outer Newton iteration method and the inner Krylov subspace iteration method. Stability analysis shows that the presented implicit weighted compact nonlinear scheme is unconditionally stable. Numerical results indicate that the implicit scheme can achieve the designed third-order accuracy in time and has a great advantage in the computation efficiency compared to the third-order explicit total variation diminishing Runge–Kutta weighted essentially non-oscillatory scheme. In addition, the implicit scheme can capture discontinuities and shock waves with high resolution and can solve Burgers’ equations with all kinds of Reynolds numbers.

High-order bicompact schemes for the quasilinear multidimensional diffusion equation

High-order bicompact schemes are designed for the quasilinear multidimensional diffusion equation. They are constructed using the method of lines, the finite volume method, the bicubic (or tricubic) Hermite interpolation. Implicit-explicit and diagonally-implicit Runge-Kutta methods are applied to time integration. The resulting schemes are stable for any ratios between grid steps, are conservative, have approximation of the fourth order in space and the third order in time. To implement the implicit-explicit schemes, an iteration method is proposed based on the approximate factorization of their multidimensional difference operators. This method is modified to implement the schemes with diagonally-implicit Runge-Kutta time stepping. High-order grid convergence and implementation efficiency of the new bicompact schemes are demonstrated on numerical examples.

Comparison of implicit time-discretization schemes for hybridized discontinuous Galerkin methods

The present study is focused on the application of two families of implicit time-integration schemes for general time-dependent balance laws of convection-diffusion-reaction type discretized by a hybridized discontinuous Galerkin method in space, namely backward differentiation formulas (BDF) and diagonally implicit Runge-Kutta (DIRK) methods. Special attention is devoted to embedded DIRK methods, which allow the incorporation of time step size adaptation algorithms in order to keep the computational effort as low as possible. The properties of the numerical solution, such as its order of convergence, are investigated by means of suitably chosen test cases for a linear convection-diffusion-reaction equation and the nonlinear system of Navier-Stokes equations. For problems considered in this work, the DIRK methods prove to be superior to high-order BDF methods in terms of both stability and accuracy.

High-order solution of Generalized Burgers–Fisher Equation using compact finite difference and DIRK methods

The main goal of this paper is to developed a high-order and accurate method for the solution of one-dimensional of generalized Burgers-Fisher with Numman boundary conditions. We combined between a fourth-order compact finite difference scheme for spatial part with diagonal implicit Runge Kutta scheme in temporal part. In addition, we discretized boundary points by using a compact finite difference scheme in terms of fourth order accuracy. This combine leads to ordinary differential equation which will take full advantage of method of line (MOL). Some numerical experiments presented to show that the combination give an accurate and reliable for solving the generalized Burgers-Fisher problems.

Development and verification of PWR core transient coupling calculation software

In PWR three-dimensional transient coupling calculation software CORCA-K, the nodal Green's function method and diagonal implicit Runge Kutta method are used to solve the spatiotemporal neutron dynamic diffusion equation, and the single-phase closed channel model and one-dimensional cylindrical heat conduction transient model are used to calculate the coolant temperature and fuel temperature. The LMW, NEACRP and PWR MOX/UO2 benchmarks and FangJiaShan (FJS) nuclear power plant (NPP) transient control rod move cases are used to verify the CORCA-K. The effects of burnup, fuel effective temperature and ejection rate on the control rod ejection process of PWR are analyzed. The conclusions are as follows: (1) core relative power and fuel Doppler temperature are in good agreement with the results of benchmark and ADPRES, and the deviation between with the reference results is within 3.0% in LMW and NEACRP benchmarks; 2) the variation trend of FJS NPP core transient parameters is consistent with the results of SMART and ADPRES. And the core relative power is in better agreement with the SMART when weighting coefficient is 0.7. Compared with SMART, the maximum deviation is −5.08% in the rod ejection condition and while −5.09% in the control rod complex movement condition.

Convergence Analysis of Hybrid High-Order Methods for the Wave Equation

We prove error estimates for the wave equation semi-discretized in space by the hybrid high-order (HHO) method. These estimates lead to optimal convergence rates for smooth solutions. We consider first the second-order formulation in time, for which we establish \({\varvec{H}}^1\) and \({\varvec{L}}^2\)-error estimates, and then the first-order formulation, for which we establish \({\varvec{H}}^1\)-error estimates. For both formulations, the space semi-discrete HHO scheme has close links with hybridizable discontinuous Galerkin schemes from the literature. Numerical experiments using either the Newmark scheme or diagonally-implicit Runge–Kutta schemes for the time discretization illustrate the theoretical findings and show that the proposed numerical schemes can be used to simulate accurately the propagation of elastic waves in heterogeneous media.

Higher‐Order Time Integration for Deformable Solids

AbstractVisually appealing and vivid simulations of deformable solids represent an important aspect of physically based computer animation. For the temporal discretization, it is customary in computer animation to use first‐order accurate integration methods, such as Backward Euler, due to their simplicity and robustness. Although there is notable research on second‐order methods, their use is not widespread. Many of these well‐known methods have significant drawbacks such as severe numerical damping or scene‐dependent time step restrictions to ensure stability. In this paper, we discuss the most relevant requirements on such methods in computer animation and motivate the interest beyond first‐order accuracy. Keeping these requirements in mind, we investigate several promising methods from the families of diagonally implicit Runge‐Kutta (DIRK) and Rosenbrock methods which currently do not appear to have considerable popularity in this field. We show that the usage of such methods improves the visual quality of physical animations. In addition, we demonstrate that they allow distinctly more control over damping at lower computational cost than classical methods. As part of our theoretical contribution, we review aspects of simulations that are often considered more intricate with higher‐order methods, such as contact handling. To this end, we derive an implicit linearized contact model based on a predictor‐corrector approach that leads to consistent behavior with higher‐order integrators as predictors. Our contact model is well suited for the simulation of stiff, nonlinear materials with the integration methods presented in this paper and more common methods such as Backward Euler alike.

A dynamically load-balanced parallel p-adaptive implicit high-order flux reconstruction method for under-resolved turbulence simulation

We present a dynamically load-balanced parallel p-adaptive implicit high-order flux reconstruction method for under-resolved turbulence simulation. The high-order explicit first stage, singly diagonal implicit Runge–Kutta (ESDIRK) method is employed to circumvent the restriction on the time step size. The pseudo transient continuation is coupled with the matrix-free restarted generalized minimal residual (GMRES) method to solve the nonlinear equations at each stage, except the first one, of ESDIRK. We use the spectral decay smoothness indicator as the refinement/coarsening indicator for p-adaptation. A dynamic load balancing technique is developed with the aid of the open-source library ParMETIS. The trivial cost, compared to implicit time stepping, of mesh repartitioning and data redistribution enables us to conduct p-adaptation and load balancing every time step. An isentropic vortex propagation case is employed to study the impact of element weights used in mesh repartitioning on parallel efficiency. We apply the p-adaptive solver for implicit large eddy simulation (ILES) of the transitional flows over a cylinder when Reynolds number (Re) is 3900 and the SD7003 wing when Re is 60000. Numerical experiments demonstrate that a significant reduction in the run time (up to 70%) and total number of solution points (up to 76%) can be achieved with p-adaptation.